# In this script we compute 5000 one-half year ahead simulations for the yield
# curve and the behaviour of our 3-bond immunized portfolio under these 
# simulated yields. We do it for Vasicek model.

# First, we run the previous script so that we get all the parameters wcomputed there.
source("S6.Fitting a model to the yield curve - CIR and Vasicek models.R")

# Parameters a, b, sigma, r(0), gamma, etc. are all in memory now.

######################## CIR Simulation #######################

set.seed(7)

# Half-year ahead
dt <- 0.5

# Nr. of simulations
n <- 5000

# We will use the affine form of the term structure:
# P(t,T) = A(t,T) * exp(B(t,T) * r(t))

# We compute A(t,T) and B(t,T) just like in script 6, but for a shift of `dt` in time:
affine.form.CIR <- tibble(`T-t` = 0:11 + dt) %>% # ...
  
# `A(t,T)` = ((2*gamma*exp((a+gamma)*`T-t`/2))/((a+gamma)*(exp(gamma*`T-t`)-1)+2*gamma))^(2*a*b/sigma^2),
# `B(t,T)` = (2*(exp(gamma * `T-t`) - 1))/((gamma + a) * (exp(gamma * `T-t`) - 1 ) + 2 * gamma)
  
affine.form.CIR

# We need to make n simulations the price of €1 for each time: P(t,t+0.5), P(t, t+1.5), ...
# First, we simulate r(t) n times
sim.r <- tibble(sim = 1:n, r_0 = r_0) %>% 
  mutate(
    # Drift term a(b-r)dt - 'drift'
    
    # Uniform(0,1) - 'unif'
    
    # dz = Normal(0,1): F^(-1)(uniform), where F is the cumulative normal function.
    # In R, this F^(-1) for the normal is qnorm().
    
    # Random term of CIR - 'rand.term'
    
    # New r(t) = r(t + dt) - 'new.r'
    
  )
sim.r

# Now we go through all the time periods and simulate the price of €1
list.sim <- vector(mode = "list", length = nrow(affine.form.CIR))

# We are going to make a cycle on the time periods (the 15 time periods considered)
for(j in seq_along(list.sim)){

  `T-t` <- affine.form.CIR %>% slice(j) %>% pull(`T-t`)
  A <- affine.form.CIR %>% slice(j) %>% pull(`A(t,T)`)
  B <- affine.form.CIR %>% slice(j) %>% pull(`B(t,T)`)
  list.sim[[j]] <- sim.r %>%
    mutate(`T-t` = `T-t`, `P(t,T)` = A * exp (-B * new.r))
  
  # Print each iteration number in the console
  cat(str_c(j," "))
  
}

list.sim[[3]] # just checking what is in tier 3 of my list (Simulations for time 3)

# Make a single table out of all 15 I have calculated and are stored in the list
sim.data.CIR <- bind_rows(list.sim)
sim.data.CIR
save(sim.data.CIR, file = "Examples/sim.data.CIR.Rdata")

# Get the Liabilities and assets cash flows (from our immunized portfolio of
# script 5)
portfolio.yield.CIR
liab.asset.cf <- immuniz.data %>%
  select(liab.cf,assets.cf,surplus) %>% 
  mutate(`T-t` = 0:11 + dt)
liab.asset.cf


sim.portfolio.CIR <- sim.data.CIR %>%
  # Add the Liabilities and assets cash flows to our simulations
  left_join(liab.asset.cf, by = join_by(`T-t`)) %>% 
  # Compute PV of each cash flow for Liabilities, assets and surplus
  #  'PV.liab.CIR', 'PV.assets.CIR', 'PV.surplus.CIR'
  mutate(
    #...
    
    
    
  ) %>% 
  # Total PV of liabilities, assets and surplus for each simulation
  group_by(sim) %>% 
  summarize(
    #...
    )

sim.portfolio.CIR
save(sim.portfolio.CIR, file = "Examples/sim.portfolio.CIR.Rdata")

# Summary of the simulation: mean
mean.surplus.CIR <- sim.portfolio.CIR %>% summarize(across(starts_with("PV"),mean))
mean.surplus.CIR

# Summary of the simulation: percentiles
q <- seq(from = 0.05, to = 0.95, by = 0.05)
quant.surplus.CIR <- sim.portfolio.CIR %>%
  # reframe() is pretty much like summary(), but for cases where the summary has
  # several lines instead of just one.
  reframe(
    quantile = q,
    percentile = scales::percent(q),
    PV.liab.CIR = quantile(PV.liab.CIR, q),
    PV.assets.CIR = quantile(PV.assets.CIR, q),
    PV.surplus.CIR = quantile(PV.surplus.CIR,q)
  )
quant.surplus.CIR

# More concise way for the last step
quant.surplus.CIR <- sim.portfolio.CIR %>%
  # reframe() is pretty much like summary(), but for cases where the summary has
  # several lines instead of just one.
  reframe(
    quantile = q,
    percentile = scales::percent(q),
    across(starts_with("PV"), ~quantile(.,q))
  )
quant.surplus.CIR

# Graphical representations: PV of surplus...
surplus.CIR.plot <- quant.surplus.CIR %>% 
  ggplot(aes(x = quantile, y = PV.surplus.CIR)) +
  geom_hline(yintercept = mean.surplus.CIR$PV.surplus.CIR, color = "darkorange") +
  geom_point(color = "darkred") +
  geom_line(color = "darkred") +
  scale_x_continuous(labels = scales::percent) +
  theme_bw()
save(surplus.CIR.plot, file = "Examples/surplus.CIR.plot.Rdata")

plotly::ggplotly(surplus.CIR.plot)

# ...and PV of projected assets and liabilities cashflows.
p <- quant.surplus.CIR %>% 
  select(!PV.surplus.CIR) %>% 
  pivot_longer(cols = starts_with("PV")) %>% 
  ggplot(aes(x = quantile, y = value, color = name)) +
  geom_point() +
  geom_line() +
  # scale_color_brewer(palette = "Set1") +
  theme_bw()
plotly::ggplotly(p)

# VaR for the PV of surplus, alpha = 99%
sim.portfolio.CIR <- sim.portfolio.CIR %>%
  # Notice that my loss variable is -PV.surplus.CIR, not PV.surplus.CIR
  # (see definition of Loss Random Variable in the slides)
  mutate(X = -PV.surplus.CIR, LossX = X - mean(X))

VaR <- sim.portfolio.CIR %>%
  summarize(VaR = quantile(LossX, probs = 0.99)) %>% 
  pull(VaR)
VaR


# TailVaR alpha = 99%

TailVaR <- sim.portfolio.CIR %>%
  filter(LossX > VaR) %>% 
  summarize(mean = mean(LossX)) %>% 
  pull(mean)

TailVaR